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ABSTRACT 

We present axisymmetric hydrodynamical simulations of the long-term accretion of a rotating gamma-ray 
burst progenitor star, a "collapsar," onto the central compact object, which we take to be a black hole. The 
simulations were carried out with the adaptive mesh refinement code FLASH in two spatial dimensions and 
with an explicit shear viscosity. The evolution of the central accretion rate exhibits phas es reminiscent of 
the long GRB 7-ray and X-ray light curve, which lends support to the proposal by Kuma r et al. ( 2008alb I 
that the luminosity is modulated by the central accretion rate. In the first "prompt" phase characterized by 
an approximately constant accretion rate, the black hole acquires most of its final mass through supersonic 
quasiradial accretion occurring at a steady rate of ~ 0.2 Mq s~^ . After a few tens of seconds, an accretion shock 
sweeps outward through the star. The formation and outward expansion of the accretion shock is accompanied 
with a sudden and rapid power-law decline in the central accretion rate M ex which resembles the Lx cx 
decline observed in the X-ray light curves. The collapsed, shock-heated stellar envelope settles into a thick, 
low-mass equatoria l disk embedded within a massive, pressure-supported atmosphere, similar to the picture 



proposed by Begel man et al.| ( [20Q8] ) for "quasistars." After a few hundred seconds, the inflow of low-angular- 
momentum material in the axial funnel reverses into an outflow from the surface of the thick disk. Meanwhile, 
the rapid decline of the accretion rate slows down, or even settles a in steady state with M ~ 5 x 10~^ Mq s~^ 
which resembles the "plateau" phase in the X-ray light curve. While the duration of the "prompt" phase 
depends on the resolution in our simulations, we provide an analytical model taking into account neutrino 
losses that estimates the duration to be ~ 20 s. The model suggests that the steep decline in GRB X-ray light 
curves is triggered by the circularization of the infalling stellar envelope at radii where the virial temperature is 
below 10^^ K, such that neutrino cooling shuts off and an outward expansion of the accretion shock becomes 
imminent; GRBs with longer prompt 7-ray emission have more slowly rotating envelopes. 

Subject headings: accretion, accretion disks — black hole physics — gamma rays: bursts — stars: winds, 
outflows — supernovae: general 



1. introduction 

Observations of long gamma-ray bursts (GRBs) carried out 
with the NASA Swift satellite have shown that the 7-ray 
prompt emission ceases after about a minute in the observer 
frame, corresponding to tens of seconds in the rest frame of 
the progenitor star. The 7-ray light curve, converted to a fidu- 
cial X-ray spectral band, smoothly joins the X-ray light curve, 
which declines rapidly, (r^ or faster) lasting for about 80 to 



curve, which was summarized by ' Zhang et al.| ( 2006[ ), reflects 
a modulation in the rate of central accretion of a rotating pro- 
genitor star onto a black hole or a neutron sta r, as in the collap- 
sar model of GRBs ( Woosley 1993; M acFadyen & Woosley| 
T999| [MacFadyen et al. 2001 ; Woosle y & Bloom|2006| ).^^ 



do not attempt to explore the implications of the potential 
presence of a magnetosphere, as in the magnetar model for 
GRBs ( e.g., Duncan & Thompson' 1992 Wheele r et al.|2000t 
Zhang & Meszaros 2001 ; Thompson et al. 2004 ; 'Komissarov 
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al.|2006| ). The rapid decline is often followed by a phase, from 
about ~ 10^ to 10"^ s, during which the X-ray flux is roughly 
constant or declines more slowly with time. The X-ray light 
curves of some GRBs exhibit "flares" where the flux increases 
suddenly by a factor of < 10^ and drops precipitously, with 
the rise and decline associated with the flare occurring on a 
time scale much shorter than the age of the burst (see, e.g.. 
Burrows et al. 2005 ; Falcone et al. 2006). Following about 
10^-1 0"^ s7a more rap id decline of the luminosity resumes 
(see, e.g., [Zhang et al.|2006| and references there in), and oc 
casionall y steepens further at ^ 10"^- 10^ s (e.g., Vaughan et 
|al.|2006| ). 

The goal of the present work is to utilize two-dimens ional 
hydrodyna mic simulations to test the hypothesis (Ku mar et 
al.|2008a|b| ) that this characteristic structure of the X-ray light 



tempt to gain insight in the origin of the steady 7-ray lumi- 
nosity (the prompt phase which we will refer to as "Phase 
0"), the rapid declin e in the X -ray light curve (Phase I in the 
nomenclature of Zha ng et al.|2 006 ), and phase of quasi-steady 
luminosity or slow decline (Phase II). We will briefly attempt 
to extrapolate the results of our simulations to the subsequent 
st eeper decline phases (P hases III and IV). 
[Kumar et al.j ( [2008a|b| ) obtained the key features of the 7- 
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ray and X-ray light curve by estimating central accretion rate 
resulting from the free (i.e., ballistic) infall of a rotating pro- 
genitor star. In this picture, the material that has sufficient ini- 
tial angular momentum to circularize outside of the innermost 
stable circular orbit (ISCO) of the black hole, forms a disk in 
the equ atorial plane, and s ubsequently accretes via disk accre- 
tion ( (Narayan et al.|200T] ). If the luminosity is then assumed 
to be proportional to the central accretion rate, and if the dis- 
tance of the 7-ray or X-ray emitting region from the center of 
the star is assumed to be approximately independent of time 
on time scales 10-10^ s, an accretion model directly trans- 
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lates into a synthetic Hg ht curve that can be co mpared with 
an observed Hght curve. |Kumar et al.l ( [2Q08a|b| ) have shown 
that with the simplest accretion model involving ballis tic in- 
fall onto the midpl ane (assumption al so made by Janiu k &| 
[Proga 2008 and Cannizzo & Gehrels ^09 ) and subsequent 
disk accretion, the mapping of the mass accretion history onto 
the light curve provides a powerful insight into the stratifica- 
tion and angular momentum str ucture of the p rogenitor star. 

In their ballistic infall model, [Kumar et aL| ( 2008a b) were 
not able to discriminate between models in which the quasi- 
steady activity in Phase II arose from disk accretion, or from 
late-time accretion from an extended stellar envelope. Depar- 
tures from ballistic infall are expected if the infalling mate- 
rial passes through an accretion shock (see, e . g.,[M acFadyen 
& Woosley"1999» 'Lee & Ramirez-Ruizi2006! 'Nagataki et al.| 
2007; Lopez-Camara et al. 2009 ), or if the disk launches a 
thermal ( MacFadyen 2003a b ; Kohri et al. 2005) or magneto- 
hydrodynamic (e.g., Proga et al. 2003) outflow ("wind") that 
can interfere with the infall. The existence of the outflow is 
particularly interesting because of the potential for nucleosyn- 
thesis in the free neutron-rich o utflow launched from the in- 
ne r part of the disk (see, e.g., jPruet et al.||20 03^, 'Surman et" 
aL][^Q 6l IFujimoto et al | [2Qa71 [Nagataki et al. 2007; Maeda 



& Tominaga 2009) anJoecause of the potential that the out- 
flow can deplete the accreting stellar envelope and limit the 
envelope mass that is accreted onto the central black hole. 

During the first ~ 10^ s following the formation of the cen- 
tral black hole when the accretion rate is M ^ 1O~^M0 s~^ 
(the precise condition depends on the black hole spin and 
shear stress-to-pressure ratio a in the disk), the inner accre- 
tion disk cools by neutrino emission and nuclear photodisin- 
tegration and accretes in a radiatively efficient fashion, except 
for in the very inner, optically thick region (e.g., Popham et 
al. 1999; Narayan et al. 2001; Di M atteo et a l. 2002; Chen 
& Belobo rodov, ,2007) . Instabilities in the thin disk have 
been cited as a candidate class of mechanisms that could 



produce th e observ ed X-ray flar es (Pema et al."2 006[ |Laz 
zati & Pema 2007 ; |Lazzati et al. 2008) and could also pro 



duce detectable gravitational radiation (P iro & Pfa hl 2007). 
Our global axisymmetric models are the necessary stepping 
stone toward the substantially more computationally demand- 
ing three-dimensional simulations that will be required to pin 
down any nonaxis ymmetric instabiliti es in the accreting col- 
lapsar (see, e.g., [Rockefeller et al.|200 6). 

We employ two-dimensional unmagnetized hydrodynamic 
simulations of the collapse, circularization, and accretion of a 
stellar envelope onto a central point mass, which we assume 
to be a black hole; relativistic corrections to the gravitational 
potential are ignored in our simulations since the innermost 
grid point lies at over 20 Schwarzschild radii in the simula- 
tion extending to the smallest radius from the black hole. The 
torque and dissipation arising from the 7? - (/) component of 
the magnetic stress is emulated with a Navier-Stokes term pa- 
rameterized by an a- viscosity prescription. For comparison 
with the X-ray light curve, we measure the central accretion 
rate. We track the flow of mass and energy at spherical radii 
10^ cm < r < 10^^ cm and interpret the results in view of the 
existing knowledge on radiatively-inefficient accretion flows. 
We observe an outflow and measure the rate at which the ac- 
creting stellar envelope is lost to the outflow. The mechanics 
of post-core-collapse accretion and outflows is key to estimat- 
ing the final mass of the black hole and the nucleo synthetic 
composition of the ejected matter (e.g., [Zhang et al.[[2008| 



and references therein). The method that we develop here 
can in future be utilized to estimate the masses of the black 
holes resulting from the collapse of massive, initially metal- 
poor "Population III" stars as well as from the collapse of the 
even more massive, hypothetical "supermassive stars," in the 
presence of rotation. 

In this work we do not simulate the neutrino-cooled disk, 
and in the simulations simply impose that the mass that 
crosses the innermost cylindrical radius of our simulation, 
^min = (0.5-2) X 10^ cm, is instantaneously incorporated in- 
side the black hole and does not provide any further ener- 
getic feedback while at radii R < R^m- This very rough 
assumption is bound to fail in general; it is most compat- 
ible with the regime in which the transition from efficient 
to inefficient cooling occurs at > ^min- Since the tran- 
sitional radius for efficient neutrino cooling recedes inward 
with the inc reasing stress-to-press ure ratio a for a given ac- 
cretion rate ( [Chen & Beloborodov|2007 ), the assumption that 
cooling is efficient within ^10^ cm is valid for a < 0.01. 
We also ignore nuclear photodisintegration when temperature 
rises above ~ (5- 10) x 10^ K; in reality, the photodissoci- 
ation allows for some heating via the capture b y free nucle- 
on s of th e neutrinos emitted in the inner disk ( [Nagataki"et| 
[aL][2007| ), which we do not model. However, we do incor- 
porate neutrino cooling in a simple analytical model for the 
evolution of the accretion shock at the radii that we do not 
resolve, ^ 5 x 10^ cm. In combination with the simulations, 
the model provides a theory for the duration of the prompt 
e mission phase observed in th e 7-rays. 

Cannizzo & Gehrels] p009| ) speculate that a cool, thin disk 
may form at large radii {R ^ 10^^ cm) at the onset of Phase 
II, and attribute the structure of the X-ray light curve to the 
long-term evolution and slow central accretion of this ex- 
tended disk. We will see that the formation of the extended 
thin disk cannot be taken for granted due to the presence of 
a massive pressure-supported convective atmosphere around 
the inner disk. 

This work is organized as follows. In Section [2j we discuss 
our numerical algorithm. In Section [3] we present the results 
our simulations. In Section |4j we present an analytical model 
for the neutrino-cooled central accretion that we do not re- 
solve in the simulations, and provide a theory for the duration 
of the prompt accretion phase and the triggering of the steep 
decline of the X-ray light curve. We also attempt to extrap- 
olate the evolution of the accretion rate beyond the duration 
of the simulations. Finally, in Section [5] we summarize our 
conclusions. 

2. NUMERICAL ALGORITHM 

The simulations were carried out with the piecewise- 
parabolic solver in the a daptive-mesh-refinement code 
FLASH ( [FryxeU et aLl|2000| ), version 2.5, in two spatial di- 
mensions using~c}Tin3ricarcoordinates (7?,z). FLASH does 
not support angular momentum adve ction and viscous trans- 
port in this regime. In Section 2.1 we describe our imple- 
mentation of angular momentum transport. In Section 2.2 we 
discuss our initial model and boundary conditions. In Section 



2.3 we provide a test of angular momentum conservation. 



2.1. Angular Momentum Transport 

The specific angular momentum £ = 7?V0, where is the 
azimuthal velocity, was treated as a mass scalar quantity that 



LINDNER ET AL. 



3 



was transported according to (see, e.g., |Pringle|1981| ) 



d(p£) 1 djRvRpi) d(v,pi 
dt dR dz 
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where is a shear viscosity to be specified below. Equation 
(l]), combined with the equation of continuity, is equivalent to 
tne azimuthal axisymmetric Navier- Stokes equation 
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(3) 



is the R-cj) component of the shear tensor. In FLASH, the 
calculation of the first three terms in equation ([T]) is carried 
out through the mass scalar advection capability; the fourth, 
parabolic term is included explicitly in the calculation of the 
radial p^-flux in the code for the diffusion of mass scalars. 

The energy dissipated through shear viscosity was ac- 
counted for by including the specific heating rate (see, e.g., 
[Landau & LiMitz|T959] ) 



R 



d_ 
dR \R^ 



(4) 



Since we do not simulate the magnetic field of the fluid, we 
utilize a local definition of the shear viscosity to emulate the 
magnetic stress arising from the intrinsically no nlocal magne- 
torotational instability (MRI; |Balbus & Hawley|1998 and ref- 
erences therein). It should be kept in mind, however, that the 
effects of MRI are in some respects very different from those 
of the viscous stress. For example, the thick disk surrounding 
our collapsar black hole is convective; in unmagnetized accre- 
tion flows convection transports angular momentum inward, 
toward the center of rotat ion ( Ryu & Goo dman] 1992[ [Stone & 



Balbus|19"96{ |Igumenshchev et al.|2000| ), whereas in magne 



tized flo ws, convection can also transport angular m omentum 
outward ([B albus & Hawley 2002 ; Igumenshchev 2002; Igu-| 
[menshchev et al. 2003 ; Christodoulou et al. 2003 ). Thus our 
results must be interpreted with caution. 

Our definition of the local viscous stress emulating the 
MRI must be valid under rotationally supported, pressure sup- 
ported, and freely falling conditions. [Thompson et al.|p005| ) 
suggest that since the wavenumber of the fastest growing MRI 
mode, which is given by the dispersion relation vaA: ~ Q 
where va is the Alfven velocity and ft = V(f)/R is the angu- 
lar velocity, should in the saturated quasi- state state be about 
the gas pressure scale height, k oc //~\ the Maxwell pv\ and 
viscous up ft stresses (up to factors in \dlnQ/dlnR\ that we 
neglect) can be equated if the viscosity is given by 



where a is a dimensionless parameter, 
height is defined locally, 

H=\V\nP\-\ 



(5) 

If the pressure scale 



(6) 



the viscosity defined in equation ([5]) suffers from divergences 
at pressure extrema. To alleviate this problem, we define a 



second viscosity according to the Shakura & Sunyaev| ( |1973 1 
prescription 



p 



(7) 



Shakura- Sunyaev viscosity overestimates the magnetic stress 
in stratified hydrostatic atmospheres. We thus set the viscosity 
in equations ([T]) and ^ to equal the harmonic mean of the 
above two viscosities 



2 ^MRI ^SS 
^MRI + ^SS 



(8) 



Our choice for the stress-to-total pressure ratio is a = 0.01, 
consistent with the ratio of the time-averaged stress to 
the time-averaged total pressure in th e stratified, radiation - 
dominated disks in the simulations ofl Hirose et al.| ( |2'009| ). 
Hirose et al. found, however, that the fluctuations in the stress 
and the pressure (total or fluid) are not temporally coincident; 
this underscores the limitations of our assumed direct propor- 
tionality of the viscous stress with the total pressure. In the 
limits umri ^ ^ss or umri <C z^ss, the effective value of the 
viscosity parameter implied by equation ^ is twice the nom- 
inal value, aeff ~ 0.02. 

Because FLASH employs an explicit method for the dif- 
fusion of mass scalars, numerical stability of the above vis- 
cous transport prescription places a stringent upper limit on 
the time step 

AR^ 

Ar<^, (9) 

where AR is the grid resolution. For a ^ 0.01, the viscous 
time step in our simulations becomes prohibitively shorter 
than the Courant time step. In our test simulations with a 7- 
law equation of state (EOS), we find that while not implying 
an outright instability, a choice of At that saturates the limit 
in equation results in weak stationary staggered perturba- 
tions in the fluid variables. We ignore this complication and 
allow our time step to be set by the limit in equation ^ of the 
cell with the smallest viscous diffusion time across the cell. 

The centrifugal force is included in the calculation of the 
gravitational acceleration via 



^grav ■ 



GMbh. fff 



(10) 



where r = (R^ -\-z^)^^^, Mbh is the mass of the central com- 
pact object which we take to be a black hole, and <^>fl(r) is 
the gravitational potential of the mass distribution of the fluid 
within the computational grid that has been spherically aver- 
aged around the origin (7?,z) = (0,0). At the radii and densi- 
ties that we resolve in the simulations, relativistic effects are 
weak; we thus treat the gravitational potential as Newtonian. 
In Section [23] below, we present a test of the angular momen- 
tum transport code. 

2.2. Initial Model and Boundary Conditions 

The initia l model is the rotating ^ 14 Mq Wolf-Rayet 
star 16TI of |Woosley & Heger] ( |2006| ), evolved to pre-core- 
collapse from a 16 main sequence progentior.^ To prepare 
the model 16TI, Woosley & Heger assumed that the rapidly 

^ Lopez-Camara et al. ( 2009 1 carried out SPH simulati ons of neutrino- 
cooled accretion during the tirst 0.5 s of the collapse and |Morsony et al!] 
(2007) simulated the propagation of a relativistic jet using the same model 
star. 
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rotating progenitor, which is near breakup at its surface at 
r ^ 4x 10^^ cm, had low metalHcity, 0.01 Z©, on the main 
sequence and became a WR star shortly after central H deple- 
tion, which implied an unusually small amount of mass loss. 
For illustration, the specific angular momentum at the three- 
quarters mass radius was 4/4 8 x 10^'^ cm^ s~\ implying 
circularization around a 10 Mq black hole at ~ 5 x 10^ cm. 
The circularization radii of the outermost layers of the star are 
in the range 10^ - 10^^ cm. Woosley & Heger provide a radius- 
dependent angular momentum profile £(r); we endowed this 
with a dependence on the polar angle 6 = cos~^(z/ r) via 



(11) 



such that the star rotates rigidly on spherical shells. 

We placed the center of the star at the origin, = (0, 0). 
Pseudo-logarithmic gridding was achieved by capping the 
adaptive resolution at radius r with A/?, Az > ^r]r, where 
we choose rj = 0.05; this prevents use of excessive resolu- 
tion far from the center of the star. Beyond the outer edge 
of the star at rstar = 4 x 10^^ cm we placed a cold (10"^ K) low- 
density stellar-wind like medium with density pro file p(r) = 



Woosley 



3 X 10"^ (r/rstar ) g cm-^ Since the model 16TI of 
|& Heger|p006| ) was not constructed to be in hydrostatic equi 
librium in the presence or rotation, we ignored rotation and 
set the initial density distribution to be spherically symmetric 
at the beginning of the simulations. This is a poor approxi- 
mation in the very outer layers of the star, as is evident from 
the ellipsoidal distortion that sets in at the beginning of the 
simulation. 

For the equation of state we chos e the Helmholtz E OS pro- 
vided with the FLASH distribution ( jTimmes & Swesty 2000 ), 
which contains contributions to pressure and internal energy 
from radiation, ions, electrons, positrons, and Coulomb ef- 
fects. We passively advected the abundances of seven nuclear 
species represented in the model including "^He, ^^C, ^^O, 
^^Ne, ^^Mg, ^^Si, and ^^Fe. The local nuclear composition 
was passed to the Helmholtz EOS as input. We do not sim- 
ulate nuclear reactions, nuclear photodisintegration, and neu- 
trino emission and absorption. These processes are certainly 
important in the hot inner accretion disk around the collap- 
sar black hole, but since we simulate only the outer, cooler 
disk with temperatures T < 10^^ K, the neglect of nuclear and 
neutrino processes is a reasonable approximation. 

The simulation was carried out in the annular cylindrical 
domain 7?min < ^ < ^max and -Zmax < z< Zmax. Bccausc the 
impact of time-step limitations (eq. ||9|) on the computational 
cost, and to avoid dealing with the fluid hot enough to be sus- 
ceptible to efficient neutrino cooling, the smallest inner ra- 
dius Rjrnn that wc could simulate was Rj^nn ~ 5 x 10^ cm. We 
placed the outer boundaries well outside the star T^max = ^max = 
(1 -5) X 10^^ cm. In TableM] we summarize the main parame- 
ters of our simulations, andalso present some of the key mea- 
surements, defined in Section |3j characterizing the outcome 
of each simulation. Each simulation was run for ~ 10^ hy- 
drodynamic time steps and consumed ~ 20,000 CPU hours 
on the Texas Advanced Computing Center's clusters Lonestar 
and Ranger. 

The boundary condition at T^min was unidirectional outflow 
that allowed free flow from larger to smaller radii (off the 
grid) and disallowed flow from smaller to larger radii (onto 
the grid) by imposing a reflecting boundary condition at Rj^i^ 
whenever vr was positive in the leftmost grid cell. We im- 
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Fig. 1 . — A test of angular momentum transport, showing the magnitude of 
the Eulerian time derivative (ri, red curve), advection (t2, green curve), and 
viscous transport (r^, blue curve) terms in the one-dimensional (d/dz = 0) 
version of the angular momentum conservation relation in equation The 
terms are expressed in the units of pi/t. We also show the sum |ri +r2 + r3 1 
{black curve), which should be much smaller than the largest of the three 
terms. The violation of angular momentum conservation near the left radial 
boundary is associated with the conservation-violating nature of the torque- 
free boundary condition that we have imposed. 

posed the torque-free boundary condition^ via 



_9_ f i_ 
dR \W 



= 0. 



(12) 



R-Rm 



As in other Eulerian codes, the boundary conditions in 
FLASH are imposed by assigning values to fluid variables in 
rows of "guard" cells just outside the boundary of the simu- 
lated domain. At any given value of z on the computational 
grid, let 7?i/2 denote the leftmost cell within the simulated do- 
main, and let Rg where g = (-| , - | , - | , - |) be the four guard 
cells to the left of 7? 1/2 such that the grid separation corre- 
sponds to = 1 . The torque-free boundary condition, if as- 
sumed to apply for R < R^mn^ implies i^jR^^- iij^jR^i^- We 
fixed the guard-cell velocity perpendicular to the left verti- 
cal boundary to vr^^ = -\vR^i/2\Ri/2/^g^ which, with the as- 
sumption of uniform density pg = pi/2, ensures mass conti- 
nuity in the guard cells and the vanishing of the mass flux 
across R = Rj^i^ if Vr i/2 > 0. All other fluid variables X were 
simply copied into the guard cells, Xg =Xi/2, and were sub- 
sequently rendered thermodynamically consistent. This sim- 
ple prescription approximates free inflow (toward smaller R) 
across Rmin, but of course, the guard cell values violate energy 
and momentum conservation atR < Rmin- 

The mass of the black hole Mbh was initialized with the 
stellar mass initially located outside the grid, at 7? < R^m- The 
black hole mass was evolved by summing mass flux crossing 
the boundary Sit R = Ru 



dMBH 

dt 



-IttR, 



/^max 
VRiRmin , Z)p(Rrmn , Z)dz. 
Zmax 



(13) 



2.3. Test of the Code 

To test our implementation of angular momentum transport, 
we performed a one-dimensional (d/dz = 0) simulation of an 



^ A motivation of the torque-free boundary condition can be found in |Zim-| 
[merman et"n!]j2005i . 
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TABLE 1 

Summary of Simulation Parameters and Key Measurements 



Run Number 


Rrnin (cm) ^ 


(/?,z)max(cm)b 


A(/?,zWn (cm)'^ 


MBH,init(M0)d 


^max (s) ^ 


a' 


fdecl (S) § 


MBH,decl(M0)h 


(JlnM/JlnOdeci ' 


1 


5x 10^ 


5 X IQii 


1.9 X 10^ 


0.51 


102 


0.01 


37 


7.35 


-2.8 


2 






6.1 X 10^ 


1.26 


103 


0.01 


47 


8.97 


-2.7 j 


3 


2x 10^ 


5 X IQii 


7.6 X 10^ 


2.05 


2x 10^ 


0.01 


52 


10.44 


-2.3 k 



^ The minimum cylindrical radius. 

^ The maximum cylindrical radius and absolute vertical latitude. 
^ The minimum resolution element size. 
^ Initial black hole mass. 

^ Duration of the simulation. 

^ The viscous stress parameter (see Section [Zl) . 

s Time of the beginning of the steep decline ot the accretion rate. 

^ Black hole mass at the beginning of the steep decline of the accretion rate. 

^ Logarithmic slope of the decline of the accretion rate. 

j For 50 s < f < 500 s. 

k For 52 s < r < 200 s. 

initially uniform temperature and density fluid with the 7 = | 
equation of state orbiting in a Keplerian potential. We use 
a uniform radial grid with T^max = 100 Rj^m and grid spacing 
AR = 0.06 Rmin- The initial temperature was chosen such that 
the sound speed was about 16% of the Keplerian velocity at 
the inner radial boundary, and 1.6 times the Keplerian veloc- 
ity at the outer radial boundary. The time step was limited by 
< ^AR^/v, which is a factor of two more stringent than 
the stability condition in equation ([9]). We found that reducing 
the time step to a half of the value required for stability sub- 
stantially reduces, but does not entirely eliminate, the noise in 
the error estimator that we are about to discuss. We evaluate 
the nonzero terms in equation ([T]) directly from the numerical 
data. Let the ri, r2, and denote the first, second, and fourth 
term in equation ([T]) 

d(pi) 

r\ = 



T2: 



T3: 



dt ' 

1 d(RvRp£) 
R dR ' 

RdR 



R^vp 



dR 



(14) 



Correct angular momentum transport requires 

|ri+r2 + r3| <Cmax(|ri|,|r2|,|r3|). (15) 

In Figure[T] we plot |ri | , |r2 1 , |r3 1 , and |ri +r2 + r3 1 over the 
entire range of radii after ~ 1 , 000 Keplerian orbital periods at 
Rrmn^ which corrcsponds to ~ 1 orbital period at T^max- Radial 
derivatives were computed by 3 -point Lagrangian interpola- 
tion with the IDL routine DERI V. In computing the derivative 
at inner boundary, we included guard cells in manner equiva- 
lent with the boundary condition prescription used in our 2D 
simulations, as described in Section |2.2| above. Some viola- 
tion of the transport equation ([T]) is expected at the two left- 
most grid cells at 7? ^ (1.03, 1.09) x R^^ because the torque- 
free boundary condition does not conserve angular momen- 
tum. Apart from the leftmost cells, the angular momentum 
is conserved at the 10% level or better at all radii. The spa- 
tial derivative in the viscous transport term (r^) is partially 
responsible for the noise evident at large radii. The noise at 
small radii seems to be correlated with the viscous time step 
limiter, which suggests that it is related to the explicit nature 
of our viscous diffusion scheme. 

3. RESULTS 



Since the evolution of the central accretion rate seems to 
fall into three distinct phases (see Figure [2]) that appear to cor- 
respond to the p hases identified in the GRB 7-ray and X-ray 
light curve (e.g., [Zhang et al.|2006| ), we divide our description 
of t he results of the simulations into three parts. In Section 

we describe Phase that concludes with the appearance 

of an accretion shock. We also discuss the limitations of our 
method in this regime, and set the stage f or an analytic model 
that we develop further below in Section |4.1| to take into ac- 



count the physics left out the our simulations. In Section 3.2 



we describe Phase I that is characterized by a steep, power- 
law decline of the central accretion rate and a rapid hydrody- 
na mic re adjustment of the accreting stellar envelope. In Sec- 
tion |33] we describe Phase II, in which the central accretion 
rate steadies. The corresponding "plateau" phase in GRB X- 
ray light curves eventually ends and gives way to a renewed 
steeper decline. Because of computational cost limitations, 
we do not extend our runs to ~ 10^ s, where, based on the 
observed light curves, one would expect the renewed steeper 
decline to occur, but further below, in Section 4.2 we briefly 
speculate about the long-term evolution of the light curve. 

3.1. Phase 0: Quasiradial Accretion 

At the radii R > T^min that we resolve in the simulations, the 
stellar collapse proceeds quasiradially for a number of sec- 
onds until an accretion shock appears at the innermost sim- 
ulated point at (R^z) = (7?min,0). The appearance of an ac- 
cretion shock coincides with the emergence of a rotationally 
supported flow in the zones at the smallest radii. In Figure 
[3j left panel, we show the density distribution, temperature 
contours, and velocity field during the quasiradial accretion 
phase, at ^ = 37 s, just prior to the formation of the accre- 
tion shock, in the simulation with T^min = 5 x 10^ cm (Run 1, 
see Table [T]). The maximum density and temperature in this 
snapshot are 2.6 x 10^ g cm~^ and T = 7.5 x 10^ K. Since 
the simulation does not allow for nuclear photodisintegration, 
the temperature in the innermost cells in the simulation is an 
overestimate. 

The existence of a quasiradial accretion phase and the late 
formation of an accretion shock are clearly artifacts of the 
choice not to simulate the innermost 5x10^ cm from the cen- 
tral axis. This innermost resolved radius is still ~ 100 times 
larger than the gravitational radius of the nascent black hole. 
Simulations that resolved the innermost radii at or near the 
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Fig. 2. — (a) Stellar mass that remains in the simulation as a function of 
time in Run 1 (blue, dot-dashed line), Run 2 (red, dashed-line), and Run 3 
(black, solid-line). The drop at r ~ 400s in Run 2 is an artifact of fluid escape 
from the box through the boundaries at = Rmax and z = ±Zmax. (b) Mass of 
the black hole, (c) The rate at which fluid mass accretes across the boundary 
at R = /?niin and is added to the mass of the black hole. The sharp drop at 
f ?^ 37-52 s coincides with the formation of the accretion shock and the onset 
of convection in the rotationally and hydrostatically supported fluid. The 
flattening at t ~ 500 s in Run 2 coincides with the cessation of the accretion 
of low angular momentum fluid through the axial funnel. The power-law 
accretion rate decline in Run 3 exhibits a shallowing of the logarithmic slope 
at t ~ 200 s. 



sarov||2QQ8[ |Komissarov & BarkQv||2QQ9} |Lopez-Camara et 
al.|2009| ), put were run much shorter than ours, saw accretion 



shock formation much earher, during the first second of the 
collapse. Some of the material falling quasiradially during 
the initial phase has enough angular momentum to circularize 
at radii that we do not resolve, but that are still larger than 
the ISCO. Indeed, our accretion shock forms earlier in runs 
with smaller T^min (see Table [T] and Figure |2]), consistent with 
the observation that circularization triggers shock formation. 
Therefore, in general, the accretion shock forms when the or- 
bital pericenter of the material crossing the equatorial plane 
becomes larger than the ISCO, or the innermost resolved ra- 
dius T^min in the s imulations in which the ISCO is unresolved. 

In Section [4?T| we will analyze differences between the in- 
ner accretion flow in our adiabatic simulations and that in the 
realistic GRBs progenitors, and suggest that the steep decline 
of the accretion rate in realistic GRB progenitors is triggered 
by the onset of circularization of the infalling stellar material 
at radii where the post- accretion- shock temperature is too low 
to allow for efficient cooling by neutrino emission. We will 
conclude that the decline seems to be associated with the on- 
set of outward expansion of the accretion shock. The outward 
expansion is distinct from and could occur much later than the 



first occurrence of the shock. In Section 4.1 we will present 
a crude analytic model in which we estimate that the trigger- 
ing of the steep decline should occur at ^deci ~ 20 s in stars 
with density and angular momentum stratification as in 16TI. 
This estimate is somewhat shorter than the shortest interval 
tdeci ^ 37 s observed in the highest-resolution simulation. Run 
1. 

3.2. Phase I: Funnel and Thick Disk Accretion 

At ^ ~ 37 - 52 s where the shortest time scale corresponds 
to the simulations that resolve the smallest radii, an accretion 
shock forms along the equator near the inner boundary (Fig- 
ure |3]) and travels outward with a velocity ~ 5 x 10^ cm s~^ 
The shocked fluid at polar angles \0-7t/2\ <75° is rotation- 
ally supported. Figure |4j left panel, shows that the isodensity 
contours of this rotationally supported fluid are roughly cir- 
cular; the vertical and the cylindrically radial pressure scale 
heights are comparable. The shocked fluid is turbulent and 
apparently convective (in two spatial dimensions, long-lived 
vortices form in the shock; the persistence of the vortices is 
an artifact of the assumed axisymmetry). The maximum tem- 
perature in the post shock fluid is 9.3 x 10^ K; this is a tem- 
perature at which the nuclear and neutrino physics that we 
ignore is marginally important; our neglect of photodisinte- 
gration cooling implies an overestimate of the temperature in 
the inner thick disk, r < 10^ cm. The shocked fluid in the 
\6-7r/2\ ^75° cone around the vertical axis continues to in- 
fall supersonically. 

Figure [5] shows the quantity p/p^, which is related to the 
specific entropy, where 7 = (Jin;?/ Jin p)^=const at ^ = 70 s and 
^ = 100 s in the highest resolution run. Run 1. Entropy appears 
to be generated throughout the thick disk. The high-entropy 
fluid exhibits a flow morphology suggestive of a "disk wind." 
The strongest outflow tracked by the highest entropy fluid is 
along the interface of the turbulent thick equatorial disk and 
the supersonic axial inflow. Prior to the cessation of the ax- 
ial inflow, the wind streamlines do not terminate at infinity, 
but rather bend back toward the equatorial plane, suggesting 
a closed meridional circulation pat tern that transports th e en- 
ergy generated in the thick disk. Ohsuga et al. p005| ) and 
[Lee & Ramirez-Ruiz ( 2006) have previously observed such a 
large-scale circulation pattern in their simulations. The high 
entropy fluid appears to accumulate at the interface of the 
thick, rotationally- supported disk and the pressure-supported 
atmosphere, and to mix convectively in the atmosphere. 

The appearance of the accretion shock is accompanied by a 
sudden rapid power-law decline of the central accretion rate. 
The times of shock formation and the onset of decline in dif- 
ferent simulations are provided in Table [T] Evolution of the 
black hole mass, the residual stellar fluia mass, and the ac- 
cretion rate, is shown in Figure |2] The steep decline of the 
accretion rate resembles the rapid decline ubiquitous in the 
observed GRB X-ray light curves. In the simulations, the de- 
cline starts at ~ 37-52 s and lasts until ~ 200-500 s. The 
logarithmic derivative of the central accretion rate during the 
decline is dlnM/dlnt ^ -2.3 -(-2.8), with the steepest de- 
cline corresponding to Run 1, the simulation that resolves the 
smallest radii. 

There does not seem to be a single explanation for the steep- 
ness of the decline of the accretion rate. We have been able 
to identify three processes that seem to contribute. We focus 
on Run 2, the run with the highest central resolution that we 
have run long enough to witness the end of the decline. 
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Fig. 3. — The innermost accretion flow, composed mostly of oxygen and neon, shortly prior (left panel) and shortly following (right panel) the formation of an 
accretion shock at ~ 37 s in Run 1 with Rmm = 5 x 10^ cm. The color scale denotes the fluid density, which in this region ranges between 10^ - 10^-^ g cm~^. 
The black contours show the T = (4, 5, 6, 7) x 10^ K isotemperature contours. The arrows show the meridional component of the fluid velocity; the longest arrows 
correspond to + vj)^/^ = 5 x 10^ cm s~^ and Mach numbers ~ 10. Following shock formation, supersonic inflow resumes along the axial funnel. 



First, the Eulerian density within the thick disk decreases 
by a factor of - 50- 100 from ^ = 50 s to 500 s. The den- 
sity drop occurs concurrently with the accretion shock expan- 
sion, and may be associated with the draining of the inner 
disk into the black hole and with a simultaneous readjusti- 
ment of the pressure-supported atmosphere of the disk toward 
near-adiabatic stratification in the presence of convection or 
large scale circulation. This decline in disk density can ex- 
plain dlnM/dlnt ^ -2 but not steeper. 

Second, there is a very gradual decrease, by a factor of < 2, 
of the vertical pressure scale height of the rotationally sup- 
ported disk during the period of the steep decline. The de- 
crease can be seen in a comparison of the left panel of Figure 
|4] showing the density distribution at ^ = 100 s, with the right 
panel of the same figure, showing the density aXt = 1,000 s. 
Since in the rotationally supported flow the viscosity is pro- 
portional to the square of the scale height, the scale height de- 
crease implies a factor of < 4 decrease of the viscosity u, and 
with it also of the disk accretion rate Mdisk- Consistent with 
the disk scale height decrease, the midplane temperature of 
the disk decreases gradually and steadily. E.g., in Run 2 at the 
innermost resolved radius of 10^ cm, the temperature drops 
from 5 X 10^ K at the onset of circularization to 2 x 10^ K at 
the end of Phase 1. 

Third, there is a rapid decline of the rate at which the low 
angular momentum fluid accretes through the axial funnel. 
Funnel accretion dominates the net accretion rate immediately 



following accretion shock formation but then drops to zero at 
the end of the steep decline at ^ ~ 500 s when the funnel in- 
flow reverses into an outflow. Our simulations may overesti- 
mate the funnel accretion rate if the funnel material is addi- 
tionally heated by a narrow relativistic axial jet, presumably 
launched from the black hole magnetosphere and responsible 
for the 7-ray and X-ray emission, that we do not simulate, but 
which must pierce the funnel region (see, e.g., Ramirez-Ruiz] 
ef al. 2002; Zhang et al. 2003, 2004; Zhang & MacFadyenj 
2006, Morsony et al. 2007; Wang et al. 2008). The effectof 
the heating of the funnel fluid by the relativistic jet might be 
to further steepen the decline of the accretion rate. Further 
magnetic outflow could develop from the corrona inner accre- 
tion disk, which could shut off funnel accretion more effec- 
tively than the outflow driven thermally by the resistive (or. 



in our approximation, viscous) dissipa tion in the disk ( |Proga 



|& Begelman|[2QQ3l [Proga et al.|[2QQ3] ). These effects could 



clearly make the central accretion rate decline, which is al- 
ready rapid in our simulations, become even more rapid. 

One also expects that the accretion rate decline is accom- 
panied by an inward recession of the boundary separating 
the radiatively-efficient, neutrino-dominated accretion flow 
(NDAF) and the radiatively-inefficent, advection or convec- 



tion dominated accretion flow (RIAF) (see, e.g., |Chen & Be- 



loborodov 2007 ). The evolution of an NDAF into an RIAF at 
radii r < 10^ cm over the course of a few hundred seconds is a 
process that may further accelerate the accretion rate decline 
in real GRB progenitors. The physics of the transition from 
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Fig. 4. — The fluid density at f = 100 s (left panel, logarithmic rendering) and t = 1,000 s (right panel, linear rendering) in Run 2. At early times, the fluid 
accreting supersonically through the axial funnel traverses multiple weak standing accretion shocks before it joins the disk or passes the boundary at = /?min- 
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Fig. 5. — The quantity p/p^ which is related to the specific entropy of the fluid, where 7 = (d\np/d\np)s=comt at constant entropy, of the fluid at r = 70 s 
(left panel) and t = 100 s (right panel) in the center of the star, in Run 1 with /^min = 5x10^ cm. The high entropy fluid tracks the outflow from the disk. The 
low-entropy fluid accreting through the axial funnel traverses multiple weak standing accretion shocks before it joints the disk or passes the boundary at = /?min- 
The primary, outward propagating accretion shock is visible along the right edge of the left panel. 



LINDNER ET AL. 



9 



o 
o 



10 



12 



" 10^° 
u 



-I— J 

o 

0) 



10' 



10' 



E 




' ' 1 1 ' 




- 

\ ^ 
















- ') 












































\ v-^- 

Vi\/\ r\ \ 








V V \ 
III 1 





10' 



10 



10 



(cm) 



X 



10 



31 



.30 



(f) 

S 10 



29 



10 



28 



1 
















\ - 




1 


1 





10' 10' 
r (cm) 



Fig. 6. — Gravitational field -d^/dr {red short dashed line), the pressure 
acceleration -p~^dp/dr {green dot-dashed line), centrifugal acceleration in 
the radial direction v^^R • r/r {blue long-dashed line), and the sum of the 
gravitational, pressure, and centrifugal acceleration {black solid line) in Run 
2 with /?min = 10^ cm. The flow is rotationally supported at r < 3 x 10^ cm 
and pressure supported at r > 4 x 10^ cm. We utilized the averages of p and 
p on spherical shells of radius r, and the mass-weighted averages of <^ and 
V0 on spherical shells. Furthermore, time-averaging was carried out in the 
intervaH = 600- 1,000 s. 



NDAF to RIAF and the onset of the outward propagation of 
the accretion shock are closely linked — bot h are controlled by 
neutrino cooling. We will argue in Section 4.1 below that the 
two operating together, starting at about the same time, is the 
most likely reason for the rather steep decline of the accre- 
tion rate between several tens of seconds and several hundred 
seconds. 

The rapid decline of the central accretion rate in our sim- 
ulations is distinct in origin fro m the le s s rapid decline seen 
in the simulations of MacFadyen et al.| (|2001i MacFadyen 
et al. simulated the fallback of the stellar envelope following 
the failure of the Shockwave resulting from the core bounce 
to unbind the star. Placing their inner numerical boundary 
at rniin,MHW = 10^ cm, they found that the radial fallback rate 
through the inner boundary declines at the rate M(rinin,MHw) 
r^/^. Since the boundary was place outside the radii of the 
infalling envelope encounters the centrifugal barrier, Mac- 
Fadyen et al. did not simulate the accretion disk and thus did 
not observe the formation and outward propagation of the ac- 
cretion shock. In our simulations, the accretion shock is aided 
by the viscous energy deposition in the rotationally-supported 
disk. The post-circularization shock seems to be responsible 
for the much more rapid central accretion rate decline in our 
simulations than in those of MacFadyen et al. 

The rapid temporal central accretion rate variability evi- 
dent in Figure [2]: is an outcome of hydrodynamical inst abili- 
ties near the innermost sim ulated radius R^^^ (see, also, |Mac-| 
Fadyen & Woosley|1999[ who observed similar variability in 
Phase 0) and should not translate into any potential variabil- 
ity of the electromagnetic jet launched from radii R <C ^min- 
The nature of the variability in the inward-directed mass flux 
M should als o be affected by the fluctuations of the magnetic 
stresses (e.g., Proga & Begelman 2003 [ Proga et aL] 2003 ) and 
by the complex interplay of the processes associated with nu- 
clear reactions and neutrino transport in the accretion flow 



Fig. 7. — The rate of mass inflow {red dashed line), mass outflow {green 
dot-dashed line), and absolute net mass flow {black solid line) crossing a 
sphere of radius r centered on the black hole, in Run 2. The rates were av- 
eraged over the time interval t = 600- 1,000 s. The outflow and the inflow 
nearly cancel over a range of radii. At the radii of the rotationally-supported 
disk r < 10^ cm, there is net inflow at the rate M ~ 5 x 10~^ Mq . 



(e.g., [MacFadyen & Woo sley'"1999^, 'Proga et al.| |2Q03t 
Igataki et al. ||2007[ |Lopez-Camara et al. 2009). TEefetere; 
we caution against ascribing phenomenological significance 
to the accretion rate variability in Figure [2]:. 

3.3. Phase II: Funnel Outflow, Thick Disk Accretion 

The steep decline of the accretion rate seems to diminish, 
or even cease at ~ 200-500 s, and the accretion rate seems to 
transition to a quasi- steady regime. This behavior resembles 
the transition into the "plateau" phase, or Phase II, of the GRB 
X-ray light curve. The simulated accretion flow appears to 
settle in a quasi- steady state, characterized by an axial outflow 
and thick equatorial disk accretion. We proceed to character- 
ize the quasi- steady accretion flow. Figure [6j which shows the 
magnitudes of the various terms in the spherically-averaged 
Euler equation and the net residual acceleration implied by 
the radial Euler equation, indicates that bulk of the fluid mass 
is rotationally supported at r < 3 x 10^ cm and is pressure 
supported at r > 4 x 10^ cm. The relative contribution of 
pressure support at radii where rotational support dominates 
is still substantial, ~ 50-75%, consistent with the thick disk 
morphology with vertical scale height h/R ^ tan 30° (Figure 
[4] right panel). Thus, our post-core-collapse accretion flow 
never resembles a thin disk. The pressure-supported atmo- 
sphere is nearly isentropic, p oc . 

In Figure [7] we show the inward-directed, the outward 
directed, ana the net mass flux flowing through spherical 
shells with radii r. After ~ 500 s, the net mass flux in 
the central ~ 10^ cm is approximately independent of ra- 
dius, which reflects a quasi- steady accretion in the inner 
part of the rotationally-supported disk at the rate M^isk ^ 
5 X 10~^ Mq s~^. This disk accreting in a quasi-steady state 
contains only Mdisk ~ 0.01 Mq, which is less than 1% of 
the mass that remains bound to the black hole in the shock- 
heated, pressure- supported atmosphere atop the rotationally- 
supported disk. The outflow and the inflow nearly cancel over 
the range of radii belonging to the atmosphere, just as was 
found in the simulat ed radiatively ine fficient accretion flow 
with convection of (Abramowicz et al.| ( |2002) . The structure 
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Fig. 8. — The large-scale density distribution {left panel) and meridional Mach number Mr-^ = (v^ + vJ)^/^/cs, where Cs is the adiabatic sound speed (right 
panel) at r = 2, 000 s in the lowest-resolution simulation, Run 3. Meridional motions in the pressure- supported atmosphere that contains most of the unaccreted 
mass are subsonic, indicating that large-scale infall has ceased. The supersonic fluid has positive Bernoulli constant and is unbound (see Figure [To|. 



described by a massive convective atmosphere surroun ding a 
thick, non radiative disk resembles the "quasistar" of Begel-" 
[man et all ( 2008 ), who envisioned the limit in which the mass 
of the pressure supported envelope exceeds the mass of the 
black hole by a large factor. 

The accretion time of the inner disk during Phase II, t^cc ^ 
^disk/^disk ~ 200 s, is shorter than duration of this phase 
(we end our simulation prior to the end of Phase II), hence 
a continuous replenishment of the inner disk must operate. 
The time scale on which the entire fluid mass bound to the 
black hole (~ 2Mq) would accrete through the thick disk is 
~ 4 X 10"^ s, though of course, not all of the mass bound at the 
beginning of Phase II must ultimately accrete; a large fraction 
could become unbound and leave in an outflow. Because of 
computational limitations we do not extend the simulations 
long enough to observe the inevitable depletion of the mas- 
sive atmosphere through inner disk accretion, but in Section 
|4.2[ we speculatively extrapolate our results into that regime. 

In Figure [sj we show a large-scale (~ 10^^ cm) view of 
the density and the meridional Mach number Mr-^ = (v| + 
v^)^/^/cs, where is the adiabatic sound speed at ^ = 2,000 s 



in the lowest-resolution simulation. Run 3. Meridional mo- 
tions in the pressure-supported atmosphere are subsonic, con- 
firming that large-scale infall has ceased. Because a vast frac- 
tion of the unaccreted mass is in this atmosphere, we nei- 
ther observe nor anticipate the tendency of the inner disk 
to spread outward in t he way in which an i s olated thin disk 



would spread and how JCannizzo & Gehrels| ( [2009| ) envision. 
At the quasi-steady disk radii, the inner and outward-directed 
mass fluxes increase outward according to Min(^), ^out(^) oc 



^1.2 



which reflects the convective or circulatory na- 
ture of the flow. The pressure-supported atmosphere at radii 
5x10^ cm ^ r < 10^^ cm contains about 0.5 Mq and exhibits 
a net inflow at the rate ^ 5 x 10""^ Mq s~\ larger than in the 
inner disk; the lack of a true steady state opens the prospect 
for a late-time, high-amplitude central accretion rate variabil- 
ity. On the other hand, the outer atmosphere r > 2 x 10^^ cm 
has a net outward-directed mass flux at 

3A/f. o-l 



' 2 



containing 

the rate Mout ~ (0.5-1.5) x 1O"^M0 s"\ though most of the 
outflowing mass remains gravitationally bound to the black 
hole. 

Figure [9] shows that the entire rotationally supported region 
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Fig. 9. — The rate of energy inflow (red short dashed line), energy outflow 
(green dot-dashed line), and absolute net energy flow (black solid line) cross- 
ing a sphere of radius r centered on the black hole in Run 2. The rates were 
averaged over the time interval t = 600- 1, 000 s. Energy outflow dominates 
inflow at all radii. The black long dashed line shows the product of pressure 
and the sound speed to which the maximum energy that can be transported 
by convection is proportional. 



2 X 10^ cm < r < 5 X 10^ cm exhibits a net outward-directed 
energy flux 

0.4 



E(r) - 10^^ erg s"^ 



2 X 109 cm 



(16) 



which in the steady state disk, r < 10^ cm, impHes a mass 
conversion efficiency ofE/ {Mc^) ~ 7 x 10""^. In the innermost 
ceUs r ~ 10^ cm, however, there is a hint of an energy inflow. 

To search for the presence of unbound flows, in Figure 
[T0| we plot the BernouUi function, defined as the sum of the 
specific kinetic ener gy, enthalpy, an d potential e nergy of the 



speciric kinetic ener gy, enthalpy, an d potential e nergy oi the 
accretion flow (e.g., [Narayan & Yi|ri 994, 19 95l|Stone et al 



1999^. 'Igumenshchev & ^Abramowicz||1999t |2000[ [Blandford 
& Begelman 1999, 2004|^ 



Be 



(17) 



where ^ is the total, negative gravitational potential. The 
structure at small radii bears only a coarse-grained resem- 
bl ance to the B ernoulli function profile in the simulations 
of |Stone et al.| ( p?9 9). The entire axial funnel region, plus 
scattered domains within the equatorial rotationally supported 
and pressure supported region, as well as the high-latitude 
fluid, |z| ^ 7 X 10^^ cm, exhibit a positive Bernoulli function. 
Be > 0, which indicates the potential for an unbound flow and 
an escape to infinity. The positivity of the Bernoulli func- 
tion does not always necessitate escape, as such fluid elements 
can be buried within the massive quasi-hydrostatic envelope 
where they can interact and mix with the Be < fluid. The 
vertically outflowing fluid above and below the rotationally 
supported disk, within an angle of ~ 20° -25° from the ver- 
tical axis, however, is entirely unbound, in clear indication 
of the presence of mass loss carried by a thermally-driven 
disk wind. Indeed, the axial outflow becomes supersonic 

^ Some authors also define the Bernoulli constant b t o equal the Bernoulli 
function divided by the Keplerian velocity, b e 
[T995; . 



Be/v^ I Narayan & Yi 



1994 



as it propagates upward through the envelope (see Fig. [8]). 
Therefore, the system spontaneously develops an advection- 
dominated inflow-outflow solution (ADIOS; |Blandford &| 
[Begelman, 1999, 2004 ). 

4. DISCUSSION 

The primary limitations of the numerical model presented 
here include: (i.) the lack of simulation coverage of the hot in- 
ner disk, r < 5 X 10^ cm, where neutrino and nuclear physics 
influence the thermodynamics of the flow, (ii.) the limited 
adequacy of the Navier-Stokes viscous fluid dynamics to ap- 
proximate the dynamics of a realistic magnetized, radiation 
dominated fluid, (iii.) the lack of modeling of the axial rel- 
ativistic jet and its enveloping cocoon, and (iv.) the lack of 
coverage of the very late evolution, ^ > 10"^ s, when the GRB 
X-ray light curve exhibits single or multiple breaks with the 
tendency toward a steepening of the luminosity decline. We 
defer an exploration of the limitations (ii.) and (iii.) for future 
work, and here briefly an d sp eculatively address limitations 
(i.) and (iv.). In Section |4.1[ we crudely take into account 
the energy loss to neutrino emission in the accretion- shock- 
heated flow and estimate the time at which the central accre- 
tion rate commences to decline steeply (the transition from 
Phase to Phase I). In Section [4!2l we discuss the implications 
of mass loss for the long-term evolution of accretion rate. 

4.1. Triggering of the Steep Decline in GRBs 

The composite 7-ray and X-ray GRB light curve starts de- 
clining abruptly and steeply after an initial period of steady 
luminosity that lasts for tens of seconds. Here, we attempt 
to shed light on the transition from the quasi- steady activity 
of the central engine to the steeply declining regime. Con- 
tinuing to work within the paradigm in which the luminosity 
is proportional to the central accretion rate, we suggest that 
the steep decline is triggered by a rapid outward expansion of 
an accretion shock through the infalling material that feeds a 
convective, rotationally- supported t hick accretion disk. This 
is consistent with the conclusion of IBamiol Duran & Kumar] 

fl08 ), who ruled out mechanisms for powering the rapid de- 
le that are intrinsic to a single emitting and cooling ele- 
ment, and by elimination inferred that the central engine (e.g., 
an accreting black hole ) remains active during t he steeply de- 



clining phase (but see Genet & Granot 2008 who showed 



that the high latitude emission from a sequence of such el- 
ements, or "pulses," could fit the decline). Because in our 
adiabatic simulations we do not resolve the innermost radii 
r < 5 X 10^ cm and do not simulate the nuclear photodisin- 
tegration, neutrino emission, neutrino capture on nucleons, 
and neutrino annihilation, that occur in the hot (T > 10^^ K) 
plasma at these radii, some of the forthcoming conclusions 
will be obtained with the aid of a simple analytical model for 
the thermodynamic evolution of the inner neutrino-cooling re- 
gion. 

We have seen that in our adiabatic simulations, the accre- 
tion shock appears when the infalling material has enough an- 
gular momentum to be rotationally supported at the innermost 
resolved radius. From there on, the accretion shock expands 
outward rapidly and sweeps through the star. The shocked 
fluid is additionally heated as a result of the viscous dissi- 
pation in the thick disk. The energy produced during the 
accretion of the thick disk is advected or convected radially 
outward by the disk "wind" (Figure |5]) which distributes it 
throughout the mostly pressure-supported torus of shocked 
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Bernoulli Function (cm^ s"^) 



5x10'^ 1x10' 
Bernoulli Function (cnn^ s"^) 



Fig. 10. — Value of the Bernoulli function Be = E]<^+Ep + [j/('y- 1)] p/p, where ^'k and Ep denote, respectively, the specific kinetic and the potential energy, 
in the simulation (left panel) and in the central 10^^ cm (right panel) at r = 1,000 s in Run 2. The material with positive values of the Bernoulli function (red 
color) has enough energy to escape from the system. 



fluid. In the adiabatic regime, after the gravitational poten- 
tial becomes dominated by the black hole while the mass 
infall rate declines rapidly, it seems that an outward expan- 
sion of the hot torus bounded by the accretion shock is in- 
evitable immediately following shock formation. In the nona- 
diabatic regime, the expansion of the shock may be delayed by 
losses associated with nuclear photodisintegration and neu- 
trino emission. If indeed, in general, a sudden and rapid drop 
in the central accretion rate accompanies the shock expansion, 
then to estimate the onset of the steep decline, one must iden- 
tify the instance at which the losses become inefficient and 
the accretion shock, aided by the viscous energy injection, 
can start traveling outward. 

Consider a fluid element with specific angular momentum 
£=10^^ iij cm^ s~^ ; its circularization radius around the black 
hole of mass Mbh = 10 Mi Mq is 



GMbh 



(18) 



If the material with density p = 10^ g cm~^ arrives at the 
shock from a free fall from infinity and the gravitational en- 
ergy density GM^up/rcirc is converted into the energy density 
in radiation aT^, where a is the radiation constant, the post- 



shock temperature is given by 



' shock ^ 



1/2 /7^xl/4 



a 



fl/2 1/4 



^6.4x lO'^K 



10x^^1 Ps 



^1/2 
-17 



(19) 



where we have taken the density jump across the shock to 
be ~ 7, appropriate if the fluid velocity is Newtonian and 
the post-shock fluid is radiation-dominated. The latter in 



particular is a crude approximation; Che n & Beloborodov 
( [2007 |) show that the pressure due to baryons and electrons 
and positrons can be comparable to and larger than the radia- 
ti on pressure in th e disk m idplane . 

|Nagataki et al.| ( |2QQ7| ) find that almost all of the energy 
emitted by neutrinos m the hot, rotationally- supported torus, 
comes from pair capture on free nucleons (t he Urea process),"^ 
for which the approximate cooling rate is ( |Qian & Woosley 



^ In the calculation of Lopez-Camara et al. ( 2009), the cooling due to neu- 
trino emission from pair annihilation dominates at early times t = 0.2 s and 
radii/? <2x 10^ cm on the equatorial plane; at later times, the cooling from 
pair annihilation is comparable to the cooling from pair capture. 
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19961 [Popham et al.|1999" and references therein) 

qeN = 9x 10'^ p8 7^10 ^nuc erg cm"^ s"^ , (20) 

where T = 10^^ Tio K is the plasma temperature, and Xnuc is 
the free nucleon fraction in statistical equilibrium, which is 
approximated via 



Xnuc ~ min[l, 8.7 p^^"^ T^^^ exp(-6.1/rio)]. 
While the cooling time 



^cool ■ 



(21) 



(22) 



where q is the net energy loss rate, which includes the pair 
capture term and other contributions, is shorter than the 
age of the collapse t, the accretion shock is confined near the 
black hole, the flow crossing the shock is highly supersonic, 
and a high accretion rate onto the central object is possible. 

When the cooling time exceeds the age of the collapse, 
^cooi > ^ the accretion shock expands outward. It seems that 
in the specific case of the collapse of a 16TI-model star, the 
pair capture neutrino emission indeed dominates losses until 
the cooling is no longer able to prevent outward expansion 
of the shock. When the post-shock temperature drops below 
~ 9 X 10^ K, which happens due to the decrease in the free 
nucleon abundance under the conditions of nuclear statistical 
equilibrium, the losses from nuclear photodisintegration and 
from neutrino emission from pair annihilation exceed those 
from pair capture, but they are not able to prevent shock ex- 
pansion. 

Consider an initial stellar density profile of the form p ex 
in which the gravity is dominated by the mass closer to the 
center, and let 



M(r) = Mo( - 



3-r 



(23) 



denote the pre-collapse mass profile, where Mq and tq denote 
the stellar mass and radius, respectively (we ignore the depar- 
ture of the density profile from a single power law near the 
stellar surface). The free fall time from radius r is given by 
tn{r) [r' /GM{r)Y/^. This relation can be inverted to obtain 
the radius, defined via %(rff) = t, from which the freely falling 
material is reaching the center at time t. 



rff(t) = (GMot^r^-')'^\ 



(24) 



The mass of the black hole grows in time and approximately 
equals 

MBii(t)r^M[rff(t)] 

r^G'/^-'Ml'^ rf-''^^t^^'/--'^ (25) 

Assuming that rdrc <C rff, the pre-shock density of infalling 
fluid at Tcirc at time t since the beginning of the explosion ap- 
proximately equals the mass infall rate ~ M[rff(t)\/t divided 
by the shock area ~ 47rr^-j.^ multiplied by the infall velocity 

^ (GMBH/rdrc)'/', 



4^[GMbh(0 4c]'/'^ 

^l-6/r ^6(l-3/r) ^12/^-5 



47r l{tf 



(26) 



We further assume that 

^(r)-4 

implying that 



M(r)' 



cm^ s ^ 



(27) 



(28) 



Kumar et al.|( |2Q08b| ) find that r ^ 2.5 throughout the bulk 
or the model 16TI ( 



of the star for the model 16TI of jWoosley & Heger| ( [2006 ) 
that we utilize and we adopt this value. We further find that 
T] ^2.5 consistent with the rotational profile of the model 
16TI in the range 4 < M(r) < 10 M© . We set ro = 10^^ cm, 
and Mo = 10 Mq and £o = 10^^-^ cm^ s~^ which approximate 
the mass and angular profile of the model 16TI, and substitute 
r = r{t) in equation p7| ) and substitute l[r{t)\ in equation (26) 
to obtain 

^9/10 33/10 



10^0 t 



-^^/^ gcm"^ 



(29) 



where in the last expression, t is given in seconds. 



Combining equations (09]), (20]), ([21]), (|^, and ([29]), the 
ratio of the cooling time for pair capture only, q = q^N, to the 
age of the collapse reads 



^cool 



2 X 10-6 ^22/5 



min[l, 17.5 ^^3/80 exp(-0.27 ^^/lO)] * 



(30) 



The ratio rises rapidly in time and becomes unity ^cooi/^ ~ 1 
at 

^decl^20s. (31) 

At this point, the material circularizing at rdrc is no longer 
able to cool by neutrino emission. Therefore, we expect that 
in a realistic star corresponding to the model 16TI, the ac- 
cretion rate starts to decline steeply at ^deci ^ 20 s after the 
explosion, when the mass of the black hole is ~ 9Mq and 
post-shock temperature is < 10^^ K. Since the initial mass of 
the black hole, if taken to equal the mass of the iron core, is 
^BH,init ~ 1.5 Mq, the implied average accretion rate preced- 
ing the decline, (M) = [MBH(Wl)-Mbh,init]/Wl ^ ^"^ 
is larger than the accretion rat e ^ (0.1-0.2) Mq s 



0.4 Mo s 
^ observed 



in our simulations a nd tho se of |MacFadyen & Woosley| ( [T999 ) 
and Nagata ki et al.|p007| , though it is consistent with the ac- 
cretion rate in the first 0.3 s in the simulation of 'Proga et al. 
( 2003 ) and in the first 0.4 s in the simulation of Lopez-Camara 
iet al.| ( |2009) . Our model possibly overestimates the infall rate 
as It does not take into account the initial hydrostatic pressure 
gradients that delay the collapse in the simulations and in real 
GRB progenitors. 

If this model for the triggering of the steep decline of 
the central accretion rate is correct, and if the onset of the 
steep decline of the accretion rate implies a termination of 
the observable prompt 7-gray emission, then the duration of 
the prompt emission in long GRBs produced by black hole- 
forming core collapse events should be anti-correlated with 
the angular momentum of the progenitor envelope. High an- 
gular momentum envelopes circularize at large radii where 
low virial temperatures imply post-shock adiabaticity and an 
earlier heating of the infalling envelope by the outward ex- 
panding accretion shock and disk outflows. Low angular mo- 
mentum envelopes may circularize at radii where the high 
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virial temperatures imply rapid cooling over several tens of 
seconds after the initial collapse. The cooling allows the ac- 
cretion shock to remain confined longer at radii where the 
free-fall velocity of the infalling envelope is highly supersonic 
and a high central accretion rate is possible, barring, of course, 
another process, such the elec t romangetically-driven outflow 



obs erved in, e.g.,|Proga et 



andlKomissarov & Barkov 



accretion. 



aD^OS]), [Nagataki et aL] (^007 ), 
l 2009| ), that could suppress central 



4.2. The Long-Term Evolution 



The period of quasi-steady or gradually declining luminos- 
ity in GRB X-ray light curves lasts for ~ 10^ - 10"^ s (Phase 

II) . At the end of this period, a steeper decline resumes (Phase 

III) , but with the shallower slope Lx oc r^-^ than in the steeply 
declining regime of Phase I. Occasionally, at r ~ 10^ - 10^ s, 
an even steeper decline Lx oc ta kes over (Phase IV in the 
nomenclature of |Zhang et al.|2QQ6| ). If the steepening of the 
light curve reflects an underlying decline of the central accre- 
tion rate, what process is responsible for this decline? Possi- 
bilities include a transformation of the character of the accre- 
tion flow due to an internal redistribution of material inside 
the accreting envelope and a depletion of the mass reservoir 
that feeds the central accretion. 

At the densities of the accreting, pressure-supported enve- 
lope, which are > 10~^ g cm~^ at the outer boundary of the 
simulation box at the end of each run, the radiation is effec- 
tively trapped and internal radiation transfer in the disk and 
the envelope is not important on the time scales on which X- 
ray light curve data are available. It also seems that given 
the near-hydrodynamic equilibrium state at ^ ~ 10^ s, any 
longer-term internal hydrodynamic redistribution of material 
between the disk and the envelope should be very gradual, so 
such a hydrodynamic redistribution is probably not a candi- 
date for the steepening that marks the transition from Phase 

II to Phase III or that which marks the transition from Phase 

III to Phase IV. Depletion of the reservoir consisting of the 
rotationally supported disk and the pressure supported at- 
mosphere could occur due to the accretion of the fluid into 
the black hole ("draining"), due to a hydrodynamic outflow 
launched from the surface of the thick disk and escaping 
through the axial funnel region ("venting"), and due to a 
radiatively-driven mass loss in the photosphere of the enve- 
lope ("blowoff"). 

In the absence of mass loss to unbound flows, the time scale 
on which the gravitationally bound envelope drains into the 
black hole, estimated from the bound envelope mass (Mgnv ~ 
2 M q) a nd the accretion rate (M ~ 5 x 10~^ Mq s~\ see Sec- 
tion[331) at r = 600 - 1 , 000 s in Run 2 is 



Me. 



M 



^4x10^ s. 



(32) 



This time scale is somewhat longer but within uncertainties 
consistent with the time scale of the initial steepening of the 
light curve at the transition from Phase II to Phase III. If 
the evolution of the envelope under draining is self- similar 
M (X Mqyw, one might expect an exponential decline of the ac- 
cretion rate; such a self- similarity, however, is not necessarily 
expected. 

The positive Bernoulli function of the fluid in the region 
of the axial funnel in Figure 10 and the mass influx and out- 
flux that increases with radius in Figure [7] suggest the pos- 
sibility that the dominant depletion may not be to accretion 



into the black hole, but instead to the loss exacted by the 
wind launched thermally from the surface of the thick, con- 
vective, rotationally supported disk. The peak net outflow rate 
at r ~ 3 X 10^^ cm in Run 2 is Mout ~ 10~^ Mq s~\ which im- 
plies a short depletion time scale of ^loss ~ 500 s. This time 
scale, however, is almost certainly an overestimate given that 
the high outflow rate may be a transient associated with the 
incomplete readjustment to the passage and breakout of the 
primary accretion shock. It seems evident, however, that the 
draining into the black hole and the mass loss to hydrody- 
namically and thermally driven outflows from the surface of 
the thick disk and the massive envelope can provide explana- 
tions of the termination of quasi- steady accretion marking the 
Phase II to Phase III transition. 

Convective energy transport in the massive, pressure- 
supported envelope can continue out to some critical radius 
where the convective motions become supersonic, resulting 
in shocks, or where radiation diffusion across the convec- 
tive cells thwarts the convective instability. Outside of this 
radius, energy is transported either radiatively or by non- 
convective bulk motions. Since the energy flux is a factor of 
10^-10^ above the Eddington limit, radiation pressure accel- 
erat es the fluid outward, resul ting in a supersonic wind (see , 
e.g., |Shaviv|2QQT]|Owocki era l. 2004; van Marie et al.|2008| ). 
The wind mass loss rate is limited by energy conservation 



l^wind^esc ^ if wind driving is radiative, momen- 

tum conservation, MwindVoo ^ L/c, where Vesc is the escape 
velocity from the critical radius, and v^o is the velocity of the 
wind at infinity. To our best knowledge, the mechanics of 
mass loss in this extremely super-Eddington regime have not 
been explored. There is the possibility that the atmospheric 
mass loss occurs on a time scale compatible with the final 
steepening of the GRB X-ray light curve, at the Phase III to 
Phase IV transition. Alternatively, as we have argued above, 
the Phase II to Phase III transition, and the Phase III to Phase 
IV transition, could both be caused by non-radiative losses 
(the draining into the black hole and venting in the axial fun- 
nel), but longer-term simulations are required to check this 
possibility. 

5. CONCLUSIONS 

We have conducted hydrodynamic simulations of the vis- 
cous post-core-collapse a ccretion of a rapidly ro tating ~ 
14 Mq Wolf-Ray et star of |Woosley & Heger[ ( [2006| ) onto the 
central black hole. The axially-symmetric simulations were 
carried out for up to 2,000 s and resolved the radii down to 
5 X 10^ cm where the collapsing stellar material circularizes 
around the black hole. The evolution of the central accretion 
rate in the simulations resembles the evolution of the observed 
GRB X-ray luminosity, which lends support to the hypothesis 
(Kumar et al. 2008a b ) that the X-ray luminosity is propor- 
tional to the rate with which stellar material accretes onto the 
black hole. We have identified three phases in the evolution 
of the accretion rate in our simulations, which appear to cor- 
respond to Phases (t he prompt phase). Ph ase I, and Phase II 
in the nomenclature of [Zhang et al.| ( |20()6] ). 

In the initial phase that in the simulations lasts 37-52 s, 
the accretion of low- angular-momentum material is quasira- 
dial for r > 5 X 10^ cm and occurs at quasi-constant rate of 
~ 0.2 Mq s~^ The end of this phase is marked by the forma- 
tion of an accretion shock at the smallest resolved radii. The 
shock immediately propagates radially outwards through the 
supersonically infalling stellar envelope. Simultaneously with 
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the formation and the outer movement of the accretion shock, 
the accretion rate drops suddenly and precipitously. We argue 
that the somewhat late onset of the accretion shock is an arti- 
fact of our not resolving the innermost two decades in radius 
outside the black hole's gravitational radius. 

We supplement the simulations with an analytical model of 
the innermost accretion disk not resolved in the simulations, 
and suggest that the accretion shock forms early, within a frac- 
tion of the first second of the formation of the black hole, 
as several published simulations of the innermost neutrino- 
cooled region have shown, but only starts to propagate out- 
ward after 20 s, when the material that is reaching the equa- 
torial plane has enough angular momentum to circularize at 
radii where the virial temperature is below ~ 10^^ K and the 
cooling by neutrino emission is suppressed. 

During the second phase characterized by a steep decline 
(X r^-^ of the accretion rate that lasts ~ 500 s, the accre- 
tion shock sweeps through the star, but a supersonic ac- 
cretion of the shocked fluid in the axial funnel region pro- 
ceeds unabated, at least in our simulations where the funnel 
has not been heated to high temperatures by the relativis- 
tic jet. The thick disk containing rotationally- supported and 
pressure- supported fluid is convective; a high-entropy out- 
flow from the inner, rotationally- supported region follows the 
accretion shock on its traversal through the star but remains 
bound within the star and appears to form a large-scale cir- 
culation pattern. The steepness of the accretion rate decline 
seems to be the consequence of a rapid hydrodynamic read- 
justment of the shocked, convective, and circulating stellar 
envelope. 

The steep decline of the accretion rate slows down or stalls 
after ~ 600 s, which appears to reflect the settling of a fraction 
of the stellar envelope in the state of near-hydrostatic equilib- 
rium. The inner, rotationally supported thick disk contains 
^ 1% of the mass of the unaccreted envelope and extends 
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to ~ 3 X 10^ cm. The thick disk is surrounded by a much 
more massive, pressure -supported atmosphere, which acts as 
a mass supply to the thick disk. At no point d o we find ev- 
idence f or the extended thin disk envisioned by ICannizzo & 
|Gehrelst (i20Q9 ). The fluid above and below the thick disk is 
mostly unbound and the simulations thus exhibit a form of a 
"disk wind." 

We speculate that depletion of the envelope through accre- 
tion onto the black hole or mass loss in thermal outflows or 
winds could be responsible for the renewed steepening of the 
GRB X-ray light curve after 10^-10"^ s. More speculatively, 
the additional steepening of the light curve occasionally ob- 
served after 10"^ - 10^ s could be due to a pervasive thermal or 
radiatively-driven mass loss in the outer layers of the atmo- 
sphere. 
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